Identification of key genes and signaling pathways related to Hetian sheep wool density by RNA-seq technology

Hetian sheep is a breed of sheep unique to the Hetian area of Xinjiang whose wool is used for producing blankets. Individual differences and hair follicle density are the key factors affecting wool production. Therefore, this study aimed to assess the Hetian sheep having different wool densities to statistically analyze the wool traits and hair follicle parameters. Furthermore, the transcriptome sequencing analysis was performed on the skins with different wool densities. The results showed that wool quantity and total hair follicle density of the high wool density sheep was significantly higher than low wool density sheep. The sheepskin with high wool density was found to grow more and finer wool than sheepskin with low wool density. A total of 1,452 differentially expressed genes were screened from the two sets of samples, including 754 upregulated and 698 downregulated genes. The differentially expressed genes were involved in the TGF-β/BMP and MAPK signaling pathways related to hair growth. Eleven differentially expressed genes belonging to the KAPs and KIFs might affect the fineness of the wool. The key genes, like the TNF, MAP2K2, INHBA, FST, PTPN11, MAP3K7, KIT, and BMPR1A, were found to probably affect the growth and density of the wool. The qPCR verified eight genes related to the MAPK pathway whose gene expression trends were consistent with the transcriptome sequencing results. This study furnishes valuable resources for enhancing the quality and production of wool in the Hetian sheep.


Introduction
Hetian sheep is a breed of wool-producing sheep that lives in the Hetian area of Xinjiang. It has the characteristics of heat resistance, roughage tolerance, and strong disease resistance. The Ministry of Agriculture has included Hetian sheep in the "National Protection List of Livestock and Poultry Genetic Resources." Hetian wool has long length, strong elasticity, high whiteness, and good bending resistance [1]. The Hetian carpet with a long history is made of Hetian wool. The carpet is famous globally because it has bright color, and a luxuriant feeling, and the carpet is comfortable, strong, and durable. The two types of hair follicles that support a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Experimental Animal Center of Tarim University. The wool density at 10 cm behind the scapula was measured by wool stemple method (Agricultural Industry Standard of the People's Republic of China, NY/T 1817-2009).
The 1 cm 2 of wool was collected five times for each sheep, collected in zip locks bags were designed to measure the wool traits. The wool density of each sheep was recorded as the average of the five measurements using the following method: The wool bundle sample is required to be neat and complete. Cut a 1 cm wool bundle as the measuring wool sample, carefully put the wool sample into a 30 mL ×60 mL weighing bottle, and wash it with petroleum ether with a boiling point of 90˚C~120˚C. Weigh the total mass of the wool sample to be tested (K mg). Take a small bunch of wool (not less than 300 roots) from each weighed wool sample, and measure the total number (N 1 roots) and total mass (K 1 mg) of the wool. Then, the formula for calculating the number of wool (N roots) with a skin area of 1 cm 2 is as follows: N = (K�K 1 )×N 1 .
Ten sheep with high wool density and ten sheep with low wool density were divided into two groups, HS means sheep in the high wool density group, and HP means sheep in the low wool density group. We selected three sheep each from HS group and HP group. The wool at 10 cm from the posterior edge of the scapula was shaved, and the skin was disinfected and injected with local anesthetics. Samples were collected using a skin punch with a diameter of 1 cm. The samples were placed into the RNase-free cryotubes and quickly frozen in liquid nitrogen.

Analysis of the wool traits
The wool samples of HS group and HP group were divided into two fiber types with fine wool and coarse wool. The wool fiber diameter of each sheep is measured. Randomly selecting 10 from the fine wool and randomly selecting 10 from the coarse wool. The diameter of the wool fiber was measured at a magnification of 100 times using a length measuring tool of the microscopic imaging system (DS-Fi3 Nikon Japan) using an optical microscope (Eclipse E200MV Nikon Japan). The morphology of the wool fiber was examined under a scanning electron microscope (EVO18 Zeiss Germany).

Morphological examination of the skin
The skin samples were fixed for 48 h, dehydrated in a series of gradient alcohol concentrations, treated with xylene, and placed in paraffin at 65˚C for 8 h. The skin samples were embedded with paraffin and were serially sectioned (cross-section; 6 μm) using a microtome (RM2016 Leica, Germany). The sections were subsequently placed on glass slides and incubated at 55˚C for 5 h. The paraffin sections were subjected to the hematoxylin-eosin (H&E) staining (following the conventional procedure) and all the images were captured using an optical microscope (Eclipse E200MV Nikon Japan) and a microscopic imaging system (DS-Fi3 Nikon Japan). The parameters of the hair follicle were determined for six sheep, and three skin slices were analyzed from each sheep. For each section, ten different microscopic fields were randomly selected from the cross-section with the sebaceous gland layer, each having an area of 1 mm 2 . The numbers of primary hair follicles and secondary hair follicles were counted in each field of view and each data is the average of the number of hair follicles counted on a section under ten different microscopic fields [23].

Transcriptome sequencing
TRIzol reagent (Invitrogen, CA, USA) was used according to the instructions to extract total RNA from the skin samples of the HS group (designated HS01, HS02, HS03) and HP group (designated HP01, HP02, and HP03). mRNA fragment was purified and reverse transcribed into cDNA. The adaptor was connected to the end of the double-stranded cDNA. The Illumina TruSeq RNA sample preparation kit (Illumina, San Diego, USA) was used to create a cDNA library by polymerase chain reaction (PCR). Six independent libraries were sequenced on the Illumina HiSeq 2500 instrument (ANOROAD, Beijing, China), according to the manufacturer's instructions.

Data analysis
FastQC was used to check the qality of the original RNA sequencing reads. The data was processed to remove impurities and obtain clean reads for subsequent bioinformatics analysis. The sieved out clean reads were mapped to the Ovis aries genome (Oar_v3.1). The raw count of each gene was determined by the default generalized linear model. Differential gene screening was performed using the factor of difference (fold change value) and q-value (adjusted pvalue, p-value after correction) as related indicators. The genes with |log2(fold change)|>0.585 and q<0.05 were considered to be differentially expressed genes (DEGs).
All the DEGs were mapped to the Gene Ontology (GO) database (http://www. geneontology.org/) with p � 0.05 as the threshold to define the GO functional classification with significant DEG enrichment. We used public pathway data (http://www.kegg.jp/) to conduct Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis on the DEGs with p � 0.05 as the threshold to define the pathways with significant DEG enrichment. Lianchuan Biological Cloud Platform (https://www.omicstudio.cn/index) was used for GO and KEGG visualization. Search tool for the Retrieval of Interacting Genes/Proteins STRING database (https:// www.string-db.org/) was used in the protein-protein interaction (PPI) analysis of the DEGs. Cytoscape 3.8 was used for producing the PPI network maps and the hub clustering modules of the PPI network were further analyzed [24].

Quantitative real-time PCR
The relative expression levels of eight DEGs screened from the DEG list were evaluated by qPCR to verify the RNA-seq results. The primer information is listed in S1 Table. Total RNA was extracted from the six skin tissues by TRIzol reagent (Invitrogen, CA, USA), and cDNA was synthesized by Reverse transcription kit (Transgen, Beijing, China). qPCR was performed in a BYQ6619-651598 instrument (Bioer Technology, Hangzhou, China) using TransStart Tip Green qPCR SuperMix (Transgen, Beijing, China). GAPDH as the internal reference gene. Amplification program: pre-denaturation at 94˚C for 30 s, denaturation at 94˚C for 5 s, and extension at 60˚C for 30 s (40 cycles). qPCR data were generated from six independent samples (n = 6). The 2 −ΔΔCt algorithm was used to calculate the qPCR results. GraphPad Prism 8 software was used to make histograms.

Statistical analyses
The statistical data of the two groups were compared through an independent sample t-test using the SPSS 23 software (SPSS Inc., Chicago, IL, USA). All the results have been expressed statistically as mean ± SD (standard deviation). P < 0.05 was considered to be a significant difference.

Characterization of wool traits
The wool of the Hetian sheep was observed to comprise fine wool and coarse wool. The wool with a diameter of fewer than 25 μm was defined as fine wool, while the wool with a diameter more than 25 μm was defined as coarse wool. The wool samples of the HS and HP groups were divided into two fiber types: fine wool and coarse wool (Fig 1A and 1B). When the microstructure of the fine wool and coarse wool was observed, the two kinds of wool fibers could be clearly distinguished under the scanning electron microscope ( Fig 1C). Then, the wool traits of the HS and HP groups were measured and analyzed statistically ( Table 1). The HS group demonstrated an average wool density of 1862.7 per cm 2 , while the HP group showed an average wool density of 1474.2 per cm 2 . The wool density of the HS group was found to be significantly higher than that of the HP group (p < 0.01). The ratio of the fine wool to coarse wool in the HS group was significantly higher than that in the HP group (p < 0.01). The average diameters of the fiber of coarse wool and fine wool in the HS group were 34.48 μm and 20.48 μm, respectively. On the other hand, the average diameters of the fiber of coarse wool and fine wool in the HP group were 40.96 μm and 21.4 μm, respectively. The diameters of the coarse wool fiber in the HS group were significantly higher than those in the HP group (p < 0.05). However, there was no significant difference in the diameter of the fine wool fiber.

Morphological analysis of the skin
The skin morphology of the Hetian sheep with different wool densities was observed. The H&E staining revealed the skin cross-section of the HS and HP groups (Fig 2A and 2B), with abundant hair follicles and sebaceous glands. The periphery of the primary hair follicles normally accompanied two sebaceous glands while the secondary hair follicles were found to be  smaller in diameter than that of the primary hair follicles, and were clustered in distribution. The densities of the total hair follicle, primary hair follicles, and secondary hair follicles, as well as the ratio of the secondary to primary hair follicles between the two groups, were analyzed statistically ( Table 2). The average total hair follicle density of the HS group was found to be 39.82 per mm 2 , and the average total hair follicle density of the HP group was found to be 31.95 per mm 2 . The hair follicle density of the HS group was found to be significantly higher than that of the HP group (p < 0.01). There was no difference in the density of the primary hair follicles between the two groups of skin. The average density of the secondary hair follicle of the HS and the HP groups were 33.79 per mm 2 and 25.93 per mm 2 , respectively. While the density of the secondary hair follicle of the HS group was found to be significantly higher than that of the HP group (p < 0.01). The ratio of the secondary to primary hair follicles in the HS group was 5.64, the ratio of the secondary to primary hair follicles in the HP group was 4.23.
The ratio in the HS group was significantly higher than that in the HP group (p < 0.01).

Distribution statistics of the expression of the gene
Six sample cDNA libraries were constructed to explore the molecular mechanism of the difference in wool density of Hetian sheep. Each library was sequenced and analyzed using Illumina HiSeq 2500, and the clean reads obtained were compared with the sheep reference genome (S2 Table). We calculated the correlation coefficient and studied the correlation to understand the overall similarity and difference of the gene expression patterns among all samples. Expectedly, all replicates showed high correlation in all comparisons ( Fig 3A). Subsequently, we screened the DEGs with |log2(fold change)|>0.585 and q < 0.05. We found 20,367 transcripts that were  jointly expressed by the two groups and 1,452 DEGs, including 754 upregulated and 698 downregulated genes, through the analysis of known coding genes ( Fig 3B). A further hierarchical aggregation of DEGs was carried out. The skin samples of the three sheep in the HS group had similar gene expression patterns, which can be clearly distinguished from the expression patterns of the samples in the HP group ( Fig 3C). The result indicates that the biological replicates have good reproducibility and the two groups of samples have differences in gene expression patterns.

GO analysis of DEGs
GO enrichment analysis revealed 1,452 DEGs were divided into 61 two-level categories. Biological processes accounted for 26 GO terms, molecular functions accounted for 14 GO terms, and cellular components accounted for 21 GO terms. The top 50 significantly enriched GO terms (p < 0.003) are shown in Fig 4. The functions of 82 genes belong to cellular developmental process, the functions of 45 genes belong to tissue development, the functions of 15 genes belong to the regulation of MAPK activity, and the functions of 13 genes belong to extracellular structure organization. In addition, the HS and HP groups had differences in BMP signaling pathway, bicellular tight junction assembly, extracellular matrix organization, activation of MAPK activity, negative regulation of MAPK activity, and bicellular tight junction assembly.

KEGG analysis of DEGs
The 1,452 DEGs were divided into six primary KEGG categories (Fig 5A), namely, organic systems, metabolism, environmental information processing, cellular processing, genetic information processing, and human diseases. The organic systems category includes 10 secondary categories, metabolism includes 12 secondary categories, environmental information processing includes 3 secondary categories, cellular processing includes 4 secondary categories, and genetic information processing includes 4 secondary categories. The key signaling pathways that affect wool grow usually belong to signal transduction pathway in environmental information processing, which includes 28 signaling pathways (Fig 5B). MAPK signaling pathway-fly, Sphingolipid, phospholipase D, MAPK, VEGF, Hippo, and TGF-β signaling pathway were found to be significantly enriched (p < 0.05). Table 3 shows the corresponding upregulated and down-regulated genes.

PPI network analysis
The PPI network analysis was performed for 169 DEGs that were enriched in twenty-eight signaling pathways in the environmental information processing (Fig 6). There were 123 proteins with existent interaction relationships and the size of each node was displayed according to the degree value (Fig 6A). A higher degree value indicated more biological processes to be involved. To identify other important clusters from the PPI network, a module analysis was conducted and 3 modules with the highest scores were selected. Cluster 1 contained 10 nodes and 18 interactions (Fig 6B). GO enrichment analysis suggested that genes in cluster 1 be mainly enriched in "protein serine phosphatase complex", "myosin phosphatase activity", "protein dephosphorylation". Cluster 2 contained 8 nodes and 12 interactions (Fig 6C), where GO enrichment analysis suggested that genes in cluster 2 to be mainly enriched in "sphingomyelin metabolic process", "ceramide biosynthetic process", "sphingolipid biosynthetic process". Cluster 3 contained 7 nodes and 10 interactions (Fig 6D), where the GO enrichment analysis suggested the genes in cluster 3 to be mainly enriched in "scaffold protein binding", "MAP kinase activity", "activation of MAPK activity", "positive regulation of MAP kinase activity". The 10 top nodes were identified by the degree value ( Fig 6E). The highly expressed genes in the HS samples included MAPK12, MAPK13, MAP2K2, TNF, PKN1. The highly expressed genes in the HP samples included RRAS2, PTPN11, BTK, VAV3, KIT. GO enrichment analysis suggested 10 genes to be mainly were enriched in the "activation of MAPK activity", "intracellular signal transduction", "positive regulation of kinase activity", "positive regulation of transferase activity", "positive regulation of MAP kinase activity". These identified genes might be partly related to the growth, and density of the wool (Fig 6B-6E).  (Table 4).

qPCR verification results
Eight genes related to the MAPK pathway were selected for qPCR verification to verify the reliability of the data. The internal reference gene was GAPDH. The results show that the gene expression data trends identified by relative quantification and RNA-seq are the same. The RNA-seq data is reliable, but the absolute expression levels from RNA-seq and qPCR are inconsistent (Fig 7). Six genes (MAPK12, RASA3, MAPK7, MAPK13, HSPB1, and MAPK2) were highly expressed in the HS group, and two genes (CHUK and EIF4E) were highly expressed in the HP group.

Discussion
Wool-producing sheep breeds from Xinjiang were used as experimental animals, and the density of Hetian wool from sheep with similar body weights and ages under the same feeding conditions was quite different. According to research, the animals with slender hair fibers were found to have a higher density of the total hair follicle than the animals with coarse hair fibers, and the fiber diameter was negatively correlated with the density of the total hair follicle [25]. A similar observation was made in the Hetian sheep. On the whole, the sheep with high wool density were found to grow finer wool and the wool quality was found to be better than the sheep with low wool density. Moreover, the total hair follicle density and the ratio of the secondary to primary hair follicles of the high wool density skin were found to be higher than that of low wool density skin. Wool fiber is the external part of the hair follicle, and the quality of wool depends on the density of the hair follicle, as well as the type and size of the hair follicle [26]. From a genetic viewpoint, the main trait affecting the "quality" and "quantity" of the wool is the follicle density [27]. This study explored the difference in the wool traits of the Hetian sheep and compared the difference in the skin hair follicle density thereby providing a better analysis of the molecular determinants of these phenotypes. Hair growth has a complex mechanism regulated by a variety of endogenous and exogenous factors [28]. Understanding the genetic basis of wool phenotype assist in improving the efficiency of wool production, and in identifying the molecular determinants of wool density and phenotype differences. At present, there are few to explore the genes and signaling pathways related to wool density, similar studies have been conducted mainly in rabbits [29][30][31]. Given this, the mRNA expression profile of the Hetian sheep skins with high and low wool densities was analyzed using RNA-seq technology, and the function of DEGs was searched to explain the molecular determinants of these differences in the wool and skin traits.
The GO analysis showed a large proportion of DEGs to be significantly enriched in the cellular developmental process, tissue development, regulation of MAPK activity, extracellular structure organization, bicellular tight junction assembly, and extracellular matrix organization. These DEGs might include the potential regulators of the density of the hair follicle and wool growth in the Hetian sheep. The KEGG enrichment analysis of the DEGs suggested seven significant enrichment pathways related to environmental information processing ( Table 3). The MAPK and TGF-β signaling are known to function in regulating hair follicle morphogenesis and hair growth. The Sphingolipid, phospholipase D, VEGF, and Hippo signaling pathways might also have a role in wool growth and development.
MAPK plays an important role in cell proliferation, differentiation, growth, aging, and death [32,33]. The MAPK signaling pathway regulates hair follicle circulation, wool growth, skin and hair follicle development, epidermal cell and keratinocyte differentiation, periodic wool growth, and wool fiber quality [34][35][36][37]. In addition, the upstream genes of the MAPK signaling pathway also play an important role in the regulation of the hair cycle and hair follicle stem cell self-renewal [38,39]. In this study, most of the DEGs related to MAPK signaling are upregulated in skin with high hair density. These up-regulated genes might have a positive effect on the density and growth of wool. According to research, the TGF-β signaling maintains the epidermal and hair cycle homeostasis, regulating the proliferation, differentiation, and apoptosis of the different epithelial stem cell populations in the hair follicles [40]. Eight DEGs were found to be enriched in the BMP signaling pathway. BMPs serve as the regulator of vertebrate development and are known to have an important role in the morphogenesis and regeneration of hair follicles [41][42][43]. In addition, BMP signaling pathway signals can also control the hair follicle cycle by regulating the proliferation and differentiation of hair matrix precursor cells; it also has a regulatory effect on the size of hair follicles [44]. BMP signaling delays the development of hair follicles and may control the distance between follicles during embryogenesis [45]. Therefore, 15 DEGs in the TGF-β/BMP signaling pathway might affect the density and growth of the wool. Most genes in the TGF-β/BMP signaling pathway were found to express higher in the sheep having low wool density. The PPI network analysis was performed on the DEGs related to environmental information processing. The 10 top nodes were identified, the up-regulated genes in the HS samples included MAPK12, MAPK13, MAP2K2, TNF, PKN1, the up-regulated genes in the HP samples included RRAS2, PTPN11, BTK, VAV3, KIT. TNF participates in hair follicle circulation and delayed hair follicle entry into catagen phase [46,47]. TNF plays an important role in wound-induced hair anagen growth and hair follicle neogenesis. In epidermal stem cells, TNF causes AKT phosphorylation and β-catenin activation [48]. Mutations in the MAP2K2 gene can cause patients to show facial deformities, skin and hair abnormalities, physical dysplasia [49]. No reports have associated MAPK12, MAPK13, PKN1, RRAS2, BTK and VAV3 to hair follicle and hair. These several genes are worthy of attention, and whether they play a role in regulating wool growth and regulating wool density requires further research. The INHBA and FST are up-regulated expressed in HS skin. INHBA is a member of the TGF-β superfamily, which plays a role in reproduction and development. Homozygous mouse lacking INHBA show dysfunctional beard, palate, and tooth development [50]. FST overexpression in sheep remarkably promotes the proliferation of sheep fetal fibroblasts and human keratinocytes [51]. In summary, the increased expression of TNF, MAP2K2, INHBA, and FST may have a positive effect on wool growth and wool density.
Mutations in the PTPN11 gene can lead to the rough hair, sparse eyebrows, skin keratinization, and obvious hypertension [52]. A study on keratinocyte-specific MAP3K7-deficient mice showed that MAP3K7 can regulate the growth, differentiation, and apoptosis of keratinocytes [53]. MAP3K7 defects can remarkably decrease E15.5 fetal mouse hair follicle precursors and delay the morphogenesis of hair follicles; the loss of MAP3K7 in hair follicles that are already in the anagen phase can lead to the passive apoptosis and degeneration of hair follicles and hair follicle damage [54]. Mutations in the KIT gene can cause a rare autosomal dominant melanocyte developmental disorder, which leads to skin and hair depigmentation and decreased hair density; patient with this disorder have thin hair, auburn hair, white hair [55,56]. The destruction of the BMP receptor, BMPR1A, in mutant mouse leads to the severe impairment of the differentiation of the inner root sheath, the reduction of hair follicles in the dermis and subcutaneous tissues, and the decline of circulating epithelial cells; BMPR1A is necessary to complete tooth morphogenesis and regulate the terminal differentiation and proliferation of hair follicles after birth [57,58]. In summary, PTPN11, MAP3K7, KIT, and BMPR1A genes may have an effect on the morphogenesis and wool growth of hair follicles, and they are all highly expressed in sheep skin with low wool density.
In addition, we found 14 DEGs belonging to KAPs and KIFs, eleven of which were up-regulated in high-wool density skin and three were up-regulated in low-wool density skin. We can come to the conclusion from the results of the wool phenotype of Hetian sheep with different wool densities, the Hetian sheep with high wool density grows more and finer wool, but the wool of low wool density Hetian sheep is thicker overall. In view of these findings, these 11 genes (LOC101116626, IFFO1, KRTAP12-1, LOC106990889, LOC105604748, KRTAP12-2,  KRTAP16-1, KRT13, KRT15, KRT2, KRTAP1-3) may play a role in affecting the fineness of wool.
Finally, we verified the sequencing results by qPCR. The results showed that the sequencing results can correctly reflect the expression ability of the detected genes. The result shows that using this technology is an economical, rapid, efficient, and reliable method for the screening of DEGs in skin tissue. Analyzing gene expression at the genome level, exploring the internal connections, and studying the molecular biological mechanisms that lead to differences in wool traits are powerful tools for studying wool quality.

Conclusions
In summary, this study has highlighted significant differences in the wool quantity and total hair follicle density in the Hetian sheep with high and low wool density. The skin with high wool density tends to grow more and finer wool than the skin with low wool density. It has detected the mRNA expression profile of the Hetian sheepskin with different wool densities and has screened out 1,452 DEGs. The DEGs involved in the TGF-β/BMP and MAPK signaling pathways related to hair growth. Eleven genes belonging to the KAPs and KIFs which might affect the fineness of the wool. The key genes, such as TNF, MAP2K2, INHBA, FST, PTPN11, MAP3K7, KIT, and BMPR1A, might affect wool growth and density. This study furnished valuable resources for enhancing the wool quality and quantity of the Hetian sheep.
Supporting information S1